1. Field of the Invention
The present invention relates to an apparatus and method for computing multiple integral, and a recording medium. More particularly, to an apparatus and method for high speed multidimensional integration of given integrand functions which are defined over a multi-dimensional unbounded domain, and a recording medium in which a program for realizing the apparatus and the method for computing multiple integral.
2. Description of the Related Art
Problems relating to computation of multiple integral on a multi-dimensional unbounded domain appear in various technical fields such as risk evaluation in financial derivative calculation, evaluating interference noises in multiple-user communication systems such as CDMA (Code Division Multiple Access), and the like.
In case of numerical integration on 1- or 2-dimensional domain, the domain (interval) can be divided into fine subdomains, and an approximated integral can be obtained as an average of represented values of an integrand function on each subdomain.
If the integration domain is in higher dimension, however, the domain must be divided into the great number of subdomains. The number of the subdomains must exponentially increase with respect to the dimensionality of the integration domain to maintain the commutative efficiency. This problem is called xe2x80x9ccurse of dimensionalityxe2x80x9d.
Monte Carlo method has been proposed as one of solutions for such the problems. Monte Carlo method usually employs pseudo-random numbers (including quasi random numbers such as low-discrepancy sequences) on the finite interval which should be uniformed as possible.
In a case where integral must be defined in an unbounded domain such as infinite domain [xe2x88x92∞, +∞], however, it is impossible to normalize uniformed random numbers on a bounded domain after converting the uniformed random numbers on the finite interval into uniformed random numbers on the unbounded domain by a linear transformation. In this case, a useful non-linear change of variables is required.
The inverse function method was established by J. von Neumann and S. Ulam in 1947. According to the method, it is able to generate random numbers in accordance with an arbitrary density function from uniformed random numbers on the finite interval. In that year, they carried out computer simulation for tracing fission neutrons, and it was successful. Since the method utilizes an inverse function, it is called inverse function method. More precisely, the inverse function xcex7xe2x88x921(X) is used for generating uniformed
xcex7(X)=∫cxcfx81(x)dx(C={u|xe2x88x92∞xe2x89xa6uxe2x89xa6X})
It has been turned out that random numbers having non-uniformed distribution on a bounded domain are available, as chaos theory have been developed recently. Chaotic dynamical mapping systems are getting more significant matters for random number generation.
Recurrence formula is usually applicable to the random number generation. In this case, a rational map is often used as the recurrence formula which generates ideally chaotic sequences. The rational map there is a result of the addition theorem of an elliptic function (including a trigonometric function). The random number generation by such the method has the following advantages.
(1) Non-cyclical random number sequence generation by which the numbers are proven to be chaotic (without some exceptions).
(2) According to sensitivity to initial conditions of chaos, completely different random number sequences are available just by slightly changing the seeds (the initial values to be applied to the recurrence formula).
(3) Known analytic function acts as the density function expressing random number distribution.
Known maps which bring the above advantages are: an Ulam-von Neumann map, a cubic map, a quintic map, and the like.
Ulam-von Neumann Map: f(x)=4x (1xe2x88x92x)
Cubic Map: f(x)=x (3xe2x88x924x)2 
Quintic Map: f(x)=x(5xe2x88x9220x+16x2)2 
Regardless of the different maps above, xcfx81(x)=1/(xcfx80(x(1xe2x88x92x))1/2) represents a density function expressing distribution of a random number sequence x[i] (ixe2x89xa70) (x[0], x[1], x[2], . . . ) which results from the following recurrence formula:
x[i+1]=f(x[i]) (ixe2x89xa70)
Japanese Unexamined Patent Application KOKAI Publication No. H10-283344 by Umeno et al. discloses a technique for generating random numbers on a bounded domain with using those maps, and embodiments where the technique is applied to Monte Carlo numerical integration. The following references disclose theoretical backgrounds for the above described application.
S. M. Ulam and J. von Neumann, Bull. Math. Soc. 53 (1947) p1120.
R. L. Adler and T. J. Rivlin, Proc. Am. Math. Soc. 15 (1964) p794.
K. Umeno, Method of constructing exactly solvable chaos, Phys. Rev. E (1997) vol. 55, pp. 5280-5284.
The above conventional technique, however, includes the following problems.
Monte Carlo numerical integration on unbounded domain requires random numbers which are distributed in the unbounded domain in accordance with a density function given by a known analytic function. Therefore, obtained uniform random numbers must be converted into the random numbers distributed in the unbounded domain. To obtain random numbers distributed in accordance with a desired density function with using the aforementioned von Neuman""s inverse function method, an integral function of the density function must be obtained first, and then, an inverse function of the integral function must be obtained.
However, it is almost impossible to analytically integrate the density function. Moreover, approximation and calculation time accompanying to numerical calculation must be considered seriously when obtaining the inverse function of the integral function resulting from the density function. Since the Monte Carlo method requires the great number of random numbers, it is a hard task to calculate the integral with using the Monte Carlo method by the inverse function method.
The numerical integration in a unbounded domain is also required for solving the Black-Scholes equation for stock option pricing based on the data of stock prices. If integration on a unbounded domain is useful for solving the Black-Scholes equation.
It has been turned out, however, that density function expressing actual distribution of stock price fluctuation shows power-law behavior (R. N. Mantegna and H. E. Stanley, Nature, vol. 376, pp. 46-49, 1995). Integral of such the density function showing the power-law behavior can not be expressed by elementary functions. Therefore, it is very difficult to solve the stock price fluctuation with using the Monte Carlo integration on an unbounded domain based on known uniform random numbers.
Such the problems also appear in evaluation of communication traffic or interference noises in CDMA showing power-law behavior in transmission time distribution or in background noise distribution respectively.
The application of process to efficient computation of multiple integral on a unbounded domain has recently been attracting attention in industry fields.
It is an object of the present invention to provide an apparatus and a method for computing multiple integral, and a recording medium storing a program for realizing the above, which are suitable for high speed numerical integration of integrands which are defined over a multi-dimensional unbounded domain.
The invention for achieving the above object will now be disclosed in accordance with the principal of the present invention.
An apparatus for computing multiple integral according to a first aspect of the present invention computes multiple integral represented by a multidimensional integrand function A to be integrated with using a vector map f. The vector map f has an unbounded support, and converts an m (mxe2x89xa71)-dimensional vector having a real number components into an m-dimensional vector having real number components. A multidimensional density function p for the limiting distribution resulting from repeatedly applying the map f to an m-dimensional vector u is analytically solvable. apparatus. Components and operations of the computing apparatus will now be described with reference to FIG. 1.
A computing apparatus comprises a first storage unit, a second storage unit, a first computing unit, a second computing unit, an update unit, and an output unit.
The first storage unit stores an m-dimensional vector x=(x1, x2, . . . , xm).
The second storage unit stores a scalar value w.
The first computing unit which computes a vector xxe2x80x2=f(x)=(xxe2x80x21, xxe2x80x22, . . . , xxe2x80x2m) resulting from applying the vector map f to the vector x being stored in the first storage unit.
The second computing unit which computes a scalar value wxe2x80x2=A(x)/xcfx81(x) based on the vector x being stored in the first storage unit and the scalar value w being stored in the second storage unit.
The update unit updates the value stored in the first storage unit by storing the vector xxe2x80x2 computed by the first computing unit on the first storage unit, and updates the value in the second storage unit by adding the scalar value wxe2x80x2 computed by the second computing unit on the second storage unit.
The output unit which computes a scalar value s=w/(c+1) based on the scalar value w being stored in the second storage unit when the number of update times by the update unit becomes c (cxe2x89xa71), and outputs the scalar value s as a numerical result of the multiple integral.
The present invention is based on the following ergodic property over multidimensional unbounded domain M:
∫MA(x)dx=∫M(A(x)/xcfx81(x))xc2x7xcfx81(x)dx=(1/N)xc2x7xcexa3j=1NA(X[j])/xcfx81(X[j]) (Nxe2x86x92∞);
X[j+1]=f(X[j]) (jxe2x89xa70);
In this case, computation errors are in the order of (1/N)1/2 regardless of the number of dimensions m even at worst. Here, N is the updating time given by c+1. According to the latest study by the inventor (K. Umeno xe2x80x9cChaotic Monte Carlo Computation: a Dynamical Effect of Random-Number Generationsxe2x80x9d, Japanese Journal of Applied Physics, March 2000), it has been clarified that the computation errors can be in the order of 1/N regardless of the number of dimensions m if a condition called Super-efficiency is satisfied. Therefore, faster integration is available.
In the computing apparatus, the scalar value stored in the second storage unit first may be a result from dividing a value resulting from applying the function A to the m-dimensional vector stored in the first storage unit first by a value resulting from applying the density function p to the m-dimensional vector stored in the first storage unit first.
In the computing apparatus, the scalar value stored in the second storage unit first may be 0, and
the output unit may compute a scalar value sxe2x80x2=w/c, and outputs the scalar value sxe2x80x2 as the numerical result of the multiple integral instead of the scalar value s.
The computing apparatus may further comprise a convergence rate obtainer, a vector map selector, and an output controller.
The convergence rate obtainer may obtain convergence rate of scalar values sequentially output by the output unit while varying the number of update times c for each of plural vector maps g1, g2, . . . , gk (kxe2x89xa72) which are prepared as the vector map f.
The vector map selector may refer to the convergence rates obtained by the convergence rate obtainer, and select a vector map gh (1xe2x89xa6hxe2x89xa6k) which shows fastest convergence rate.
The output controller may control the output unit to output the scalar value with using the vector map gh as the vector map f and the number of update times cxe2x80x2 (cxe2x80x2 greater than c) instead of the number of update times c.
limiting distribution of a vector sequence
u, f(u), f(f(u)), f(f(f(u))), . . .
resulting from applying the vector map f to a predetermined m-dimensional vector u=(u1, u2, . . . , um) for equal to or greater than 0 times, may satisfy the following property of:
xcfx81(u)=Πi=0mxcfx81i(ui);
xcfx81i(ui)xcx9ccxe2x88x92i|ui|xe2x88x92(1+a) for uixe2x86x92xe2x88x92∞;
xcfx81i(ui)xcx9ccxe2x88x92i|ui|xe2x88x92(1+a) for uixe2x86x92+∞;
(a greater than 0, 1xe2x89xa6ixe2x89xa6m, cxe2x88x92i greater than 0, c+i greater than 0)
In the computing apparatus, the vector map f may be defined as
f(u)=(f1(u1), f2(us), . . . , fm(um))
by a function fi(t)=gi(di t)/di di greater than 0) which is defined in 1xe2x89xa6ixe2x89xa6m, and the map gi may be defined by any one of the following maps xcfx86j (1xe2x89xa6jxe2x89xa68) and a natural number ni (nixe2x89xa72), as follows:
gi(xcfx86j(xcex8))=xcfx86j(nixcex8);
xcfx861(xcex8)=xe2x88x92sgn(tan xcex8)/|tan xcex8|1/a;
xcfx862(xcex8)=xe2x88x92sgn(tan xcex8)xc3x97|tan xcex8|1/a;
xcfx863(xcex8)=xe2x88x92sgn(cos xcex8)/|tan xcex8|1/a;
xcfx864(xcex8)=xe2x88x92sgn(cos xcex8)xc3x97|tan xcex8|1/a;
xcfx865(xcex8)=sgn(cos xcex8)/|tan xcex8|1/a;
xcfx866(xcex8)=sgn(cos xcex8)xc3x97|tan xcex8|1/a;
xcfx867(xcex8)=sgn(sin xcex8)/|tan xcex8|1/a;
xcfx868(xcex8)=sgn(sin xcex8)xc3x97|tan xcex8|1/a;
sgn(t)=1 for t greater than 0;
sgn(t)=0 for t=0;
sgn(t)=xe2x88x921 for t less than 0
The computing apparatus may comprise a convergence rate obtainer, a positive
The convergence rate obtainer may define the map f for each of plural positive numbers q1, q2, . . . , qk (kxe2x89xa72) prepared as an invariable a, and obtain convergence rates of the scalar values sequentially output by the output unit while varying the number of update times c.
The positive number selector may refer the convergence rates obtained by the convergence rate obtainer, and select a positive number qh (1xe2x89xa6hxe2x89xa6k) which shows the fastest convergence rate.
The output controller may define the map f with using the positive number qh as the invariable a, and control the output unit to output the scalar values with using the number of update times cxe2x80x2 (cxe2x80x2 greater than c) instead of the number of update times c.
The computing apparatus may comprise a convergence rate obtainer, a map selector, and an output controller.
The convergence rate obtainer may define the map gi with using plural ones of the maps xcfx86j, and obtain convergence rates of the scalar values sequentially output by the output unit while varying the number of update times c.
The map selector may refer to the convergence rates obtained by the convergence rate obtainer, and select one of the maps xcfx86j which shows the fastest convergence rate.
The output controller may define the map g1 with using the map xcfx86j selected by the map selector, and control the output unit to output the scalar values with using the number of update times cxe2x80x2 (cxe2x80x2 greater than c) instead of the number of update times c.
The computing apparatus may comprise a convergence rate obtainer, a natural number selector, and an output controller.
The convergence rate obtainer may define the map gi relating to each of plural natural numbers p1, p2, . . . , pk (kxe2x89xa72) as the natural numbers ni, and obtain convergence rates of the scalar values sequentially output by the output unit while varying the number of update times c. convergence rate obtainer, and select a natural number ph (1xe2x89xa6hxe2x89xa6k) which shows the fastest convergence rate.
The output controller may define the natural number map gi with using the natural number ph as the natural number ni, and control the output unit to output the scalar values with using the number of update times cxe2x80x2 (cxe2x80x2 greater than c) instead of the number of update times c.
The output unit of the computing apparatus may compute the scalar value s each time the update unit carries out update, compare the latest scalar value s with the former scalar value which is computed at former update, and output the latest scalar value s if a result of the comparison satisfies a predetermined condition for terminating the computation.
A method for computing multiple integral according to a second aspect of the present invention computes multiple integral of a multidimensional integrand function A to be integrated with using a vector map f, a first storage unit which stores an m-dimensional vector x=(x1, x2, . . . xm), and a second storage unit which stores a scalar value w.
The vector map f has an unbounded support, and converts m (mxe2x89xa71)-dimensional vector having real number components into m-dimensional vector having real number components. A multidimensional density function xcfx81 for the limiting distribution resulting from repeatedly applying the map f to m-dimensional vector u is analytically solvable.
The computing method comprises the first computing step, the second computing step, the updating step, and the outputting step.
The first computing step computes a vector xxe2x80x2=f(x)=(xxe2x80x21, xxe2x80x22, . . . xxe2x80x2m) resulting from applying the vector map f to the vector x being stored in the first storage unit.
The second computing step computes a scalar value wxe2x80x2=A(x)/xcfx81(x) based on the vector x being stored in the first storage unit and the scalar value w being stored in the
The updating step updates the value stored in the first storage unit by storing the vector xxe2x80x2 computed by the first computing unit on the first storage unit, and updating the value in the second storage unit by adding the scalar value wxe2x80x2 computed by the second computing unit to a value to be stored on the second storage unit.
The outputting step computes a scalar value s=w/(c+1) based on the scalar value w being stored in the second storage unit when the number of update times by the update unit becomes c (cxe2x89xa71), and outputs the scalar value s as a numerial result of the multiple integral.
In the computing method, the scalar value stored in the second storage unit first may be a result from dividing a value resulting from applying the function A to the m-dimensional vector stored in the first storage unit first by a value resulting from applying the density function xcfx81 to the m-dimensional vector stored in the first storage unit first.
In the computing method, the scalar value stored in the second storage unit first may be 0, and
the outputting step may compute a scalar value sxe2x80x2=w/c, and output the scalar value sxe2x80x2 as the numerical result of the multiple integral instead of the scalar value s.
The computing method may further comprises the convergence rate obtaining step, and the vector map selecting step.
The convergence obtaining step may obtain convergence rate of scalar values sequentially output by the outputting step while varying the number of update times c for each of plural vector maps g1, g2, . . . , gk (kxe2x89xa72) which are prepared as the vector map f.
The vector map selecting step may refer to the convergence rates obtained by the convergence rate obtainer, and select a vector map gh (1xe2x89xa6hxe2x89xa6k) which shows fastest convergence rate.
A vector map gh may be used as the vector map f.
The outputting step may output the scalar value with using the vector map gh as the times c.
In the computing method, a multidimensional density function xcfx81 representing the limiting distribution of a vector sequence
u, f(u), f(f(u)), f(f(f(u))), . . .
resulting from applying the vector map f to a predetermined m-dimensional vector u=(u1, u2, . . . , um) for equal to or greater than 0 times, may satisfy the following property of:
xcfx81(u)=Πi=0mxcfx81i(ui);
xcfx81i(ui)xcx9ccxe2x88x92i|ui|xe2x88x92(1+a) for uixe2x86x92xe2x88x92∞;
xcfx81i(ui)xcx9ccxe2x88x92i|ui|xe2x88x92(1+a) for uixe2x86x92+∞;
(a greater than 0, 1xe2x89xa6ixe2x89xa6m, cxe2x88x92i greater than 0, c+i greater than 0)
In the computing method, the vector map f may be defined as
f(u)=(f1(u1), f2(us), . . . , fm(um))
by a function fi(t)=gi(di t)/di (di greater than 0) which is defined in 1xe2x89xa6ixe2x89xa6m, and the map gi is defined by any one of the following maps xcfx86j (1xe2x89xa6jxe2x89xa68) and a natural number ni (nixe2x89xa72), as follows:
gi(xcfx86j(xcex8))=xcfx86j(nixcex8);
xcfx861(xcex8)=xe2x88x92sgn(tan xcex8)/|tan xcex8|1/a;
xcfx862(xcex8)=xe2x88x92sgn(tan xcex8)xc3x97|tan xcex8|1/a;
xcfx863(xcex8)=xe2x88x92sgn(cos xcex8)/|tan xcex8|1/a;
xcfx864(xcex8)=xe2x88x92sgn(cos xcex8)xc3x97|tan xcex8|1/a;
xcfx865(xcex8)=sgn(cos xcex8)/|tan xcex8|1/a;
xcfx866(xcex8)=sgn(cos xcex8)xc3x97|tan xcex8|1/a;
xcfx867(xcex8)=sgn(sin xcex8)/|tan xcex8|1/a;
xcfx868(xcex8)=sgn(sin xcex8)xc3x97|tan xcex8|1/a;
sgn(t)=1 for t greater than 0;
sgn(t)=0 for t=0;
sgn(t)=xe2x88x921 for t less than 0
The computing method may further comprise the convergence rate obtaining step, the positive number selecting step, and the output controlling step.
The convergence rate obtaining step may define the map f for each of plural positive numbers q1, q2, . . . , qk (kxe2x89xa72) prepared as an invariable a, and obtain convergence rates of the scalar values sequentially output by the output unit while varying the number of update times c.
The positive number selecting step may refer to the obtained convergence rates, and select a positive number qh (1xe2x89xa6hxe2x89xa6k) which shows the fastest convergence rate.
The output controlling step may define the map f with using the positive number qh as the invariable a, and control output of the scalar values with using the number of update times cxe2x80x2 (cxe2x80x2 greater than c) instead of the number of update times c.
The computing method may further comprises the convergence rate obtaining step, the map selecting step, and the output control step.
The convergence rate obtaining step may define the map g1 with using plural ones of the maps xcfx86j, and obtain convergence rates of the scalar values sequentially output by the output step while varying the number of update times c.
The map selecting step may refer to the obtained convergence rates, and select one of the maps xcfx86j which shows the fastest convergence rate.
The output controlling step may define the map gi with using the selected map xcfx86j selected, and control output of the scalar values with using the number of update times cxe2x80x2 (cxe2x80x2 greater than c) instead of the number of update times c.
The computing method may further comprises the convergence rate obtaining step, the natural number selecting step, and the output controlling step.
The convergence rate obtaining step may define the map gi relating to each of plural natural numbers xcfx811, xcfx812, . . . , xcfx81k (kxe2x89xa72) as the natural numbers ni, and obtain convergence rates of the scalar values sequentially output by the outputting step while varying the
The natural number selecting step may refer to the obtained convergence rates, and select a natural number xcfx81h (1xe2x89xa6hxe2x89xa6k) which shows the fastest convergence rate.
The output control step may define the natural number map gi with using the natural number xcfx81h as the natural number ni, and control output of the scalar values with using the number of update times cxe2x80x2 (cxe2x80x2 greater than c) instead of the number of update times c.
In the computing method, the outputting step may compute the scalar value s each time the updating step carries out update, compares the latest scalar value s with the former scalar value which is computed at former update, and output the latest scalar value s if a result of the comparison satisfies a predetermined condition for terminating the computation.
A program for realizing the apparatus and method for computing multiple integral according to the present invention may be stored in a recording medium such as a compact disc, a floppy disk, a hard disk, a magneto-optical disk, a digital versatile disc, a magnetic tape, or a semiconductor memory.
The program stored in the recording medium according to the present invention may be executed by a data processor being equipped with a storage device, a computing unit, output device, and the like, such as a general purpose computer, a video game device, a PDA (Personal Data Assistants), and a cellular phone, to realize the above described apparatus and method for computing multiple integral.
The recording medium storing the program according to the present invention may be distributed or merchandized being separated from the data processor.